Computational Training: running fastR

Great R resources

What is “fast code”?

1. ✏️ Easy to write
   familiarity with code editor, libraries

2. 💡 Easy to understand
   structured, with consistent variable names, commented.

3. 🔧 Easy to debug
   clear naming, DRY, tests.

4. 🏎️ Easy to run
   (profiling, C++, using “optimized” code).

Variable naming

tmp <-  10
tmp1 <- tmp * 24

cases_per_hour <- 10
cases_per_day <- cases_per_hour * 24

Embrace functional programming I

Functions are first-class citizens in R

first-class-citizenship.R

f <- function(x){x^2}

lapply(1:10, f)

generator <- function(n=2){
    function(x){x^n}
}
cube <- generator(3)

list(one_function = f)

Embrace functional programming II

Rethink for, while loops; “apply” instead

Say you want to extract the \(R^2\) from three linear models with different predictors (or formulae).

formulae.R

formulae <- c(
    Sepal.Length ~ Sepal.Width,
    Sepal.Length ~ Petal.Length,
    Sepal.Length ~ Species
)

bad way

lm_results2 <- c()

for (formula in formulae) {
    fit <- lm(formula, data = iris)
    r2 <- summary(fit)$r.squared
    lm_results2 <- c(lm_results2, r2)
}

good way

extract_r2 <- function(formula) {
    fit <- lm(formula, data = iris)
    r2 <- summary(fit)$r.squared
    return(r2)
}

lm_results <- sapply(formulae, extract_r2)

side effects

exists("fit")

The parallel package

You can imagine wanting to run each of the apply/for loop iterations in parallel.

sockets (not on Windows)

library(parallel)
f <- function(i) {
    lme4::lmer(
        Petal.Width ~ . - Species + (1 | Species),
        data = iris)
}

system.time(save1 <- lapply(1:100, f))
##    user  system elapsed
##   2.048   0.019   2.084
system.time(save2 <- mclapply(1:100, f))
##    user  system elapsed
##   1.295   0.150   1.471

forking

num_cores <- detectCores()
cl <- makeCluster(num_cores)
system.time(save3 <- parLapply(cl, 1:100, f))
#    user  system elapsed 
#   0.198   0.044   1.032 
stopCluster(cl)

Introduction to data.table

• data.table is a package that extends the data.frame class.
• (Quicker) alternative to dplyr for large datasets.
• cheatsheet

all you need to know

• use the data.table called dt …
• subset it on the rows specified by i…
• and manipulate columns with j…
• grouped according to by.

library(data.table)
iris_dt <- as.data.table(iris)

iris summaries

Subsetting and Summarizing

iris_dt[Species == "setosa", mean(Sepal.Length)]

[1] 5.006

Grouping and Aggregating

iris_dt[, mean(Sepal.Length), by=Species]

      Species    V1
1:     setosa 5.006
2: versicolor 5.936
3:  virginica 6.588

iris summaries (.VARS)

Counting entries (.N)

iris_dt[Sepal.Length>5, .N, by=Species]

      Species  N
1:     setosa 22
2: versicolor 47
3:  virginica 49

Print 1st entry of each group (.SD)

iris_dt[, .SD[1], by=Species]

      Species Sepal.Length Sepal.Width Petal.Length Petal.Width
1:     setosa          5.1         3.5          1.4         0.2
2: versicolor          7.0         3.2          4.7         1.4
3:  virginica          6.3         3.3          6.0         2.5

Modifying the dt (in place)

You will need a walrus (operator):

Modifying the dt (in place)

Define a new column

iris_dt[, Species_ita := paste0(
    Species, 'ino'
), by=Species]
iris_dt[, unique(Species_ita)]

[1] "setosaino"     "versicolorino" "virginicaino" 

Modify existing column

iris_dt[, Species_ita := gsub(
    "aino","ino",
    Species_ita
), by=Species_ita]
iris_dt[, unique(Species_ita)]

[1] "setosino"      "versicolorino" "virginicino"  

Acting on multiple columns

Apply a function to multiple columns

numeric_cols <- names(iris_dt)[sapply(iris_dt, is.numeric)]
iris_dt[, lapply(
    .SD, as.integer
), .SDcols = numeric_cols] |> head(2)

   Sepal.Length Sepal.Width Petal.Length Petal.Width
1:            5           3            1           0
2:            4           3            1           0

Define multiple columns at once

new_cols <- paste0(numeric_cols, "_int")
iris_dt[, (new_cols) := lapply(
    .SD, as.integer
), .SDcols = numeric_cols]; head(iris_dt, n=2)

   Sepal.Length Sepal.Width Petal.Length Petal.Width Species Species_ita
1:          5.1         3.5          1.4         0.2  setosa    setosino
2:          4.9         3.0          1.4         0.2  setosa    setosino
   Sepal.Length_int Sepal.Width_int Petal.Length_int Petal.Width_int
1:                5               3                1               0
2:                4               3                1               0

Summary and advanced topics

Profiling

• Identify (and hopefully fix!) bottlenecks in your code.
• The profvis package is a good package to use for this purpose.

Profiling Example: Column Means

library(profvis)
library(data.table)
n <- 4e5
cols <- 150
data <- as.data.frame(x = matrix(rnorm(n * cols, mean = 5), ncol = cols))
data <- cbind(id = paste0("g", seq_len(n)), data)
dataDF <- as.data.table(data)
numeric_vars <- setdiff(names(data), "id")

profvis({
  means <- apply(data[, names(data) != "id"], 2, mean)
  means <- colMeans(data[, names(data) != "id"])
  means <- lapply(data[, names(data) != "id"], mean)
  means <- vapply(data[, names(data) != "id"], mean, numeric(1))
  means <- matrixStats::colMeans2(as.matrix(data[, names(data) != "id"]))
  means <- dataDF[, lapply(.SD, mean), .SDcols = numeric_vars]
})

Good coding practice: Reproducibility & Generalisability

Why are code reproducibility & generalisability important?

• ✅ Transparency & Verification
• 🤝 Collaboration & Longevity
• 🔍 Quicker detection of errors

Reproducibility

• Main idea: Be able to reproduce results to ensure they are valid.
• Common practices:
  • Setting the seed
  • Importing all necessary packages
  • Documenting R environment (sessionInfo())

bad way

x <- rnorm(1)
dtruncnorm(x, -5, 5, 0, 1)

still bad...

set.seed(1234)
x <- rnorm(1)
dtruncnorm(x, -5, 5, 0, 1)

there we go!

library(truncnorm)
set.seed(1234)
x <- rnorm(1)
dtruncnorm(x, -5, 5, 0, 1)

Generalisability

• Code should be “generalisable” meaning that anyone else can refer to it and use it on their own data.
• Common practices:
  • ⚙️ Write functions for operations you use frequently
  • 💬 Document your code (add comments)
  • ❓ Add basic checks (testing)
  • ✔️❌ Test your code on different data sets

Generalisability: Bad example

don't even bother...

kmeans_recip <- function(){
  for (i in c(1:4)){
    iris[, i] <- 1/iris[, i]
  }
  kmeans_res <- kmeans(iris[, c(1:4)], centers = 3)
  return(kmeans_res$cluster)
}
kmeans_recip()[1:20]

 [1] 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 1 1 1 1 1

Generalisability: OK example

• The previous example works but it’s really not generalisable…
• Obvious things that need improvement:
  • Hardcoded values (first 4 columns are numeric in iris)
  • Data set should be an input argument

it's getting better

kmeans_recip <- function(data, cont_cols){
  for (i in cont_cols){
    data[, i] <- 1/data[, i]
  }
  kmeans_res <- kmeans(data[, cont_cols], centers = 3)
  return(kmeans_res$cluster)
}
kmeans_recip(data = iris, cont_cols = c(1:4))[1:20]

 [1] 3 3 3 3 3 1 1 3 3 3 3 3 3 3 3 1 1 1 1 1

Generalisability: Can we generalise?

• Now we can also check if this works on the diamonds data set from the ggplot2 package.

seems like it works!

library(ggplot2)
kmeans_recip(data = diamonds)[1:20]

Error in ifelse(any(data[, i] == 0), stop("Division by 0 not allowed!"), : Division by 0 not allowed!

Debugging

• 🐛 Debugging is the process of finding and fixing errors in your code.
• 🤬 Generally annoying but can be made much easier with a few simple steps!

Debugging made easy

• Common practices to debug your code more easily include:
  • 📖 Read error messages carefully
  • 🖨️ Use print messages (cat(), print(), message() etc.)
  • 💾 Save objects that may be causing the error (to be able to reproduce it)
  • 📑 A bit more “advanced” options are traceback(), browser(), debug() (more details here)
  • 💻 Don’t forget: Google, Stack Overflow, ChatGPT etc. are your friends

Error Handling

• 😰 Sometimes errors cannot be avoided…
• 🎉 Good news is we can handle them!
• 💊 Error handling is the process of responding to and recovering from error conditions.

Error Handling: the tryCatch() syntax

• tryCatch() is the function to use for error handling in R.

basic tryCatch syntax

result <- tryCatch({
  expr
  },
  warning = function(w){
    warning-handler-code
    },
  error = function(e){
    error-handler-code
    },
  finally = {
    cleanup-code
  }
)

Error Handling: basic example

tryCatch example

safe_log <- function(x){
  result <- tryCatch({
    log(x) # Attempt to calculate the logarithm
  },
  warning = function(w){
    message("A warning occurred: ", w) # Handle warnings
    NULL # Return NULL if a warning occurs
  },
  error = function(e){
    message("An error occurred: ", e) # Handle the error
    NA  # Return NA if an error occurs
  },
  finally = {
    # This block executes no matter what
    message("Logarithm attempt completed.")
  })
  return(result)
}

Error Handling: basic example

tryCatch works!

print(safe_log(2024))

Logarithm attempt completed.

[1] 7.612831

print(safe_log(-8))

A warning occurred: simpleWarning in log(x): NaNs produced

Logarithm attempt completed.

NULL

print(safe_log("zero"))

An error occurred: Error in log(x): non-numeric argument to mathematical function

Logarithm attempt completed.

[1] NA